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ABSTRACT 


This  research  addresses  the  issue  of  discriminant  transportability,  that  is,  the  development 
of  discriminants  which  can  be  applied  to  events  in  a  variety  of  geographic  areas.  The  failure  of 
discriminants  to  work  for  events  in  new  areas  of  study  is  most  often  attributed  to  the  effects  of 
changes  in  the  propagation  path  (Q,  phase  blockage,  scattering,  focusing/defocusing,  etc.).  An 
additional  factor  is  the  instability  of  spectral  ratio  measurements.  The  spectra  of  Pn  and  Lg  waves 
for  small  explosions  often  exhibit  spectral  nulls  due  to  the  time-delayed  multiple  shots.  Often  these 
spectral  nulls  occur  at  slightly  different  frequencies,  contributing  to  artificial  maxima  and  minima  of 
the  Lg/Pn  spectral  ratio.  In  our  approach,  we  deconvolve  the  source  effects  separately  for  Lg  and 
Pn  before  the  spectral  ratio  is  computed.  This  is  accomplished  by  detrending  the  array  averaged 
log-spectra  and  then  computing  a  source  time  function  based  on  the  spectral  nulls.  We  use  both  the 
deconvolved  spectral  ratio  and  the  derived  source  time  function  for  event  identification.  This 
procedure  has  been  tested  on  a  number  of  Ground  Tmth  Databases  which  include  events  from 
Scandinavia,  central  Europe,  Spain,  and  the  Middle  East. 

CTBT  RELEVANCE 

Transportability  of  discriminants  is  an  important  issue  in  the  monitoring  of  a  global  CTBT 
because  the  events  which  generate  the  most  interest  wiU  likely  be  those  which  occur  in  geographic 
areas  which  have  yet  to  be  studied  with  seismic  methods.  Additionally,  if  the  scatter  in 
discriminants  can  be  reduced,  their  rehability  for  source  identification  even  in  areas  already  studied 
should  improve. 

KEY  WORDS 

Discrimination,  mining,  explosions,  signal  processing,  discriminant  transportability,  ground-truth 
databases 
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Background 


For  over  ten  years  it  has  been  known  that  there  are  observable  differences  in  the  regional 
phases  Pn,  Sn,  and  Lg  for  earthquakes  and  explosion  sources  (e.g.  Dysart  and  Pulli  (1987), 
Baumgardt  and  Ziegler  (1988),  Hedlin  et  al.  (1990)).  These  differences  were  noted  in  the 
(relatively)  high-frequency  band  of  the  data,  using  events  recorded  at  the  newly  installed 
NORESS,  ARCESS,  and  FINESA  arrays.  Specifically,  it  was  noted  that  the  Lg/Pn  and  Sn/Pn 
spectral  ratios  above  10  Hz  provided  excellent  discrimination  capability.  Research  then  focused  on 
a  variety  of  methods  to  parameterize  and  process  these  observed  differences,  including  statistical 
methods  (Shumway,  1996),  neural  network  methods  (Dysart  and  Pulli,  1990),  and  fuzzy  logic 
(Sereno  and  Wahl,  1993). 

However,  when  these  discriminants  were  applied  to  events  in  new  areas  of  interest,  their 
ability  to  identify  source  types  was  shown  to  be  inferior  when  compared  to  events  from  the  original 
source  area.  In  hindsight  this  was  not  unexpected,  given  the  area-to-area  variation  in  propagation 
effects  which  can  change  the  spectral  content  of  the  regional  phases.  Some  of  these  effects,  such  as 
anelastic  attenuation,  can  be  independently  measured  and  modeled  for  application  to  the  observed 
spectra.  Other  effects,  such  as  Lg  blockage  and  the  focusing/defocusing  of  regional  phases  due  to 
velocity  perturbations,  are  more  difficult  to  correct.  Scattering  also  contributes  to  spectral  variation 
(Dainty,  1996),  although  recent  developments  in  the  measurement  of  scattered  waves  (coda)  show 
promise  for  source  identification  (Mayeda  and  Walter,  1996). 

An  additional  consideration  is  the  actual  comparison  of  the  Pn,  Sn,  and  Lg  phase  spectra. 
Spectral  ratios  are  typically  used  for  spectral  comparison  and  these  formed  the  basis  of  our  early 
studies.  To  parameterize  the  spectral  ratios,  we  computed  the  mean  values  of  Lg/Pn  and  Sn/Pn 
over  a  variety  of  frequency  bands:  0-20  Hz  to  characterize  the  broadband  spectra  as  well  as  2-5  Hz, 
5-10  Hz,  and  10-20  Hz.  A  closer  examination  of  this  calculation  illustrates  the  difficulty  of  the 
estimation.  Fig.l  shows  the  Lg  and  Pn  spectra  for  an  explosion  at  the  Koashva  Nostoch  apatite 
mine  in  the  Apatity  area.  The  Lg  spectrum  peaks  around  2  Hz,  then  quickly  decays  by  a  factor  of 
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10,  after  which  there  is  a  smooth  decay  with  frequency.  The  Pn  spectrum  however  shows  evidence 
of  spectral  scalloping  as  is  often  seen  in  quarry  explosions.  This  is  presumably  due  to  ripple  firing 
(if  the  firing  pattern  is  delayed  by  enough  time  to  appear  in  the  pass  band  of  the  data)  or  the  total 
length  of  the  firing  pattern.  When  we  calculate  the  Lg/Pn  spectral  ratio,  these  peaks  and  troughs  of 
the  spectra  result  in  an  estimate  which  shows  simply  artifacts  of  the  spectral  ratio  and  not  a  true 
comparison  of  the  relationship  of  the  spectral  energy  in  each  phase. 


Figure  1.  Lg  and  Pn  spectra  (left)  as  well  as  Lg/Pn  spectral  ratio  (right)  for  an  explosion  at  the  Koashva  Nostoch 
apatite  mine  in  the  Apatity  area  (GTDB  event  number  2066).  The  spectral  ratio  has  been  parameterized  as  mean 
values  across  the  entire  band  as  well  as  from  2-5  Hz,  5-10  Hz,  and  10-20  Hz.  In  this  case,  the  parameterization  is 
clearly  not  representative  of  the  data. 


Approach 

In  our  approach,  we  deconvolve  the  source  effects  separately  for  Lg  and  Pn  before  the 
spectral  ratio  is  computed.  This  is  accomplished  by  first  detrending  the  array  averaged  log-spectra 
of  Lg  and  Pn,  then  computing  a  source  time  function  based  on  the  spectral  nulls.  This  source  time 
function  is  then  deconvolved  from  the  observed  spectra  in  the  frequency  domain.  The  scatter  in  the 
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resulting  spectral  ratio  parameterizations  is  then  reduced.  Additionally,  the  deconvolved  spectral 
ratios  can  be  used  in  conjunction  with  the  derived  source  time  function  as  a  combined  discriminant. 
This  procedure  has  been  tested  on  a  number  of  Ground  Truth  Databases  which  include  events  from 
Scandinavia,  Spain,  and  the  Middle  East. 

Software  and  Data  Processing 

Software  used  during  this  research  effort  has  been  written  in  the  MATLAB®  technical 
computing  environment.  The  software  consists  of  scripts  or  functions  (so-called  m-files)  which 
call  the  built-in  MATLAB®  routines  to  perform  basic  vector  and  matrix  computations.  These  scripts 
can  be  linked  to  compiled  C  or  FORTRAN  objects,  or  can  be  translated  to  ANSI  C  using  a 
supplied  compiler.  Toolboxes,  which  combine  numerous  m-files,  are  also  available  which  perform 
groups  of  functions,  such  as  signal  processing,  image  processing,  and  neural  computing.  This 
computing  environment  was  chosen  to  take  advantage  of  signal  processing  routines  already 
developed  for  other  projects,  and  because  the  binary  data  files  {mat-files)  and  ASCII  m-files  can  be 
moved  from  platform-to-platform  (UNIX,  PC,  Macintosh)  without  translation. 

The  starting  point  is  a  CMR  Version  3.0  database  {.wfdisc,  .origin,  and  .w  files;  see 
Anderson  et  al,  1990).  The  C-program  csslmatlab  reads  a  user-specified  .wfdisc  file  and  generates 
mat-files  for  each  .wfdisc  entry  with  file  names  station_channel.mat  (e.g.  ARA0_sz.mat).  Each 
mat-file  contains  variables  with  names  corresponding  to  those  in  the  CMR  Version  3.0  definition, 
as  well  as  other  comments  which  are  used  by  subsequent  plotting  routines.  At  this  point,  the 
programs  gram.m  and  zoom.m  are  used  to  display  a  single  channel  of  data  and  select  windows  for 
the  noise,  Pn,  Sn,  and  Lg  phases.  The  program  extract_arcess.m  (as  well  as  programs  for  the 
NORESS,  GERESS,  FINESA,  and  ILPA  arrays)  are  used  to  extract  the  waveform  segments  and 
produce  matrix  variables  of  the  phase  waveforms  across  each  array.  For  the  regional  data  of 
interest  to  this  project,  waveform  segments  consist  of  512-points  for  each  channel,  corresponding 
to  12.8  seconds  for  the  typical  40  samples/second  digitization.  Once  the  waveforms  have  been 
extracted,  event  files  (e.g.  elOOl.mat)  are  created  which  contain  all  of  the  relevant  information, 
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including  that  from  the  .origin  file  if  available.  Event  names  correspond  to  the  .wfdisc  prefix  name 
in  order  to  remain  consistent  with  the  Ground-Truth  Database  convention  (Grant  et.  al,  1993). 

Signal  processing  is  accomplished  with  a  variety  of  developed  m-files  which  typically 
consist  of  four  versions;  compute,  plot,  process,  and  hardcopy.  For  example,  there  is 
compute _spectra.m,  plot_spectra.m,  process_spectra.m,  and  hardcopy_spectra.m.  The  compute 
version  performs  the  actual  calculation,  the  plot  version  plots  the  results,  the  process  version 
batches  all  event  files  for  the  compute  version  and  saves  the  result  to  each  event  file,  and  the 
hardcopy  version  batches  all  plots  to  a  specified  printer,  m-files  have  been  developed  to  perform  a 
variety  of  these  signal  processing  functions,  including  spectra,  power  spectral  density,  cepstra, 
coherence,  cross-spectra,  multiple-filter  analysis,  spectral  ratios,  and  beamforming,  as  well  as 
three-component  analysis. 

Databases  Examined 

Our  earlier  studies  of  regional  seismic  event  identification  utilized  the  parameterization  of 
the  Pn,  Sn,  and  Lg  spectra  discussed  above  (mean  values  over  broad  and  narrow  bands). 
Waveform  windows  of  the  regional  phases  were  not  saved  as  part  of  the  effort,  just  the 
parameterizations  of  the  spectra.  In  order  to  address  this  new  research,  we  dearchived  all  of  the  old 
data  and  used  the  software  described  above  to  re-pick  phase  windows  and  create  event  files 
consisting  of  the  waveform  segments  of  pre-event  noise,  Pn,  Sn,  and  Lg  phases.  These  databases 
include: 

•  NORESS  events  from  1985  -  1988,  earthquakes  &  mining  explosions 

•  Ground  Truth  Database  1:  Vogtland  earthquakes  &  explosions 

•  Ground  Truth  Database  2:  Steigen  earthquake  swarm 

•  Ground  Tmth  Database  8:  Apatity  mining-induced  tremors 

•  Ground  Tmth  Database  9:  Apatity  mine  blasts 

•  Spanish  Ground  Tmth  Database:  Alhucemas  &  Puertollano  events 
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•  Galilee  Ground  Truth  Database:  earthquakes  &  mining  explosions 

Observations  and  Analysis 

Fig.  2  shows  Pn  recordings  at  the  ARCESS  array  of  six  small  explosions  at  the  Apatity 
mining  area.  In  each  case  the  array  data  have  been  steered  toward  the  source  to  illustrate  the 
correlation  of  multiple  arrivals  in  the  Pn  wave  train.  Note  that  for  Events  2022,  2023,  2024,  2039, 
and  2046,  the  Pn  wave  consists  of  2  or  3  distinct  arrivals,  whereas  for  Event  2045  the  Pn  arrival  is 
a  single  pulse.  Since  these  events  are  from  the  same  area,  this  difference  in  arrival  structure  is 
likely  a  source  effect,  resulting  from  different  firing  patterns  for  each  explosion  at  the  quarry. 

Fig.  3  shows  the  corresponding  spectra  of  Pn,  Lg,  and  pre-event  noise  for  these  same 
mining  explosions.  Note  that  for  event  2039  the  spectral  nulls  for  Pn  and  Lg  waves  occur  at  the 
same  frequencies,  but  for  event  2046  the  spectral  nulls  are  absent  for  Lg.  Fig.  4  shows  the  array 
averaged  autocorrelations  of  Pn  and  Lg.  Note  that  these  autocorrelations  indicate  a  source  time 
function  that  is  different  for  Pn  versus  Lg. 

The  objective  of  the  deconvolution  is  to  remove  the  effect  of  the  multiple  explosions  while 
retaining  the  information  from  the  individual  shot.  This  can  be  accomplished  in  a  number  of  ways. 
One  way  is  to  manually  pick  the  peaks  and  troughs  of  the  Pn  and  Lg  spectra  and  sort  the  delays. 
The  delay  pattern  is  determined  from  either  the  inverse  of  the  peak  frequency  of  the  inverse,  or  2 
times  the  trough  frequency.  A  second  way  is  to  pre-whiten  the  observed  spectra  with  a  frequency- 
squared  correction,  threshold  the  spectra,  inverse  Fourier  transform,  and  compare  the  results  with 
the  autocorrelation.  A  third  way  is  to  precompute  the  interference  spectra  for  a  wide  range  of 
possible  source  time  functions  (Fig.  5)  and  perform  a  pattern  match  between  the  models  and  the 
pre- whitened  observed  spectra.  This  is  followed  by  an  inverse  Fourier  transform  and  a  comparison 
with  the  autocorrelation. 
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GTDB  Event  2022  -  Steered  Pn 


GTDB  Event  2023  -  Steered  Pn 
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Figure  2.  Steered  Pn  waveforms  at  ARCESS  for  six  mining  explosions  at  the  Apatity  area. 
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Figure  3.  ARCESS  array  estimates  of  Pn,  Lg,  and  noise  spectra  for  six  mining  explosions  at 

Apatity. 
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GTDB  Event  2024:  Pn(-),  Lg(-)  GTDB  Event  2039:  Pn(-),  Lg(-) 


Figure  4.  ARCESS  array  estimates  of  Pn  and  Lg  autocorrelations  for  six  mining  explosions  at 

Apatity. 
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Figure  5.  Model  spectra  of  a  quarry  blast  source  consisting  of  multiple  charges. 


Fig.  6  shows  the  resulting  effect  on  the  broadband  Pn/Lg  spectral  ratios  for  100  events 
recorded  at  the  NORESS  array.  Source  type  was  determined  from  a  number  of  local  bulletin 
sources  as  explained  in  Dysait  and  Pulli  (1990).  The  upper  figure  shows  the  spectral  ratios 
calculated  without  source  deconvolution.  Note  that  the  scatter  in  this  discriminant  is  large  for 
explosions,  since  these  events  show  the  largest  variation  in  source  time  function.  Once  this  source 
time  function  is  deconvolved  from  the  spectra,  the  resulting  spectral  ratios  show  much  less  scatter. 


Incorporation  of  Derived  Source  Time  Function 

The  source  deconvolution  just  described  essentially  removes  the  information  about  the 
multiple  nature  of  the  source  in  order  to  derive  the  underlying  single  source  spectral  ratio. 
However,  the  multiple  nature  of  the  source  is  in  itself  an  indicator  of  source  type.  So  a  valid 
question  is  whether  or  not  we  can  combine  this  information  for  source  identification. 

Fig.  7  shows  the  derived  source  time  function  durations  for  the  102  NORESS  events 
analyzed  in  the  previous  section.  The  earthquakes  in  this  data  set  have  source  time  functions  that 
range  from  zero  (or  the  resolution  of  the  data)  to  0.5  seconds.  The  explosions  have  source  time 
functions  that  range  from  zero  to  1.8  seconds. 

Fig.  8  plots  these  derived  source  time  functions  versus  the  deconvolved  broadband  Pn/Lg 
spectral  ratios.  The  use  of  these  two  parameters  provides  complete  event  type  separation  for  this 
data  set. 


Extension  to  Other  Ground  Truth  Databases 

The  overall  objective  of  this  work  has  been  the  determination  of  so-called  transportable 
seismic  discriminants,  and  the  approach  has  been  to  stabilize  the  most  widely  used  seismic 
discriminant,  the  spectral  ratio.  At  this  point,  we  examine  the  use  of  this  stabilized  spectral  ratio 
along  with  the  derived  source  time  function  for  events  outside  the  area  of  initial  study.  This  has 
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o:  Earthquakes  x:  Explosions 


After  Source  Deconvolution 


Figure  6.  Broadband  Pn/Lg  spectral  ratios  at  NORESS  before  (top)  and  after  (bottom) 

deconvolution. 
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Derived  Source  Time  Function 


o:  Earthquakes  x:  Explosions 


Event  Number 


Figure  7.  Derived  source  time  function  for  earthquakes  and  explosions  recorded  at  NORESS. 


12 


Broadband  Pn/Lg 
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Figure  8.  Broadband  Pn/Lg  spectral  ratios  after  source  deconvolution  versus  derived  source  time 

function  for  NORESS  events. 
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involved  the  analysis  of  new  Ground  Truth  Databases  compiled  as  part  of  the  sponsor’s  overall 
research  program.  For  example,  the  Spanish  Dataset  includes  recordings  of  mining  explosions  and 
earthquakes  at  the  Sonseca  Array,  and  we  can  use  this  array  to  examine  whether  or  not  the 
observation  of  multiple  arrivals  in  the  Pn  wave  train  can  be  seen,  as  is  often  seen  at  the 
Scandinavian  arrays.  Fig.  9  shows  the  steered  Pn  arrivals  at  the  Sonseca  array  for  four  mining 
explosions  in  the  Puertollano  subset  of  the  database.  Athough  not  as  prominent  as  with  mining 
explosions  at  Apatity,  mulitple  arrivals  can  still  be  seen  for  these  events. 

Fig.  10  shows  the  broadband  Pn/Lg  spectral  ratio  versus  derived  source  time  function  for 
200  events  comprising  the  original  NORESS  database,  the  Apatity  mine  explosions,  the  Steigen 
aftershock  series,  the  Puertollano  events,  and  events  from  the  Galilee  earthquake  &  mining 
explosion  database.  The  earthquakes  in  this  combined  data  set  have  source  time  functions  that 
range  from  zero  (or  the  resolution  of  the  data)  to  1.1  seconds  (a  greater  range  than  for  th  hmited 
NORESS  dataset).  The  explosions  have  source  time  functions  that  range  from  zero  to  1.8 
seconds.  The  broadband  deconvolved  Pn/Lg  spectral  ratio  for  the  explosions  is  larger  than  for  the 
earthquakes,  and  the  separation  between  the  two  classes  is  nearly  complete. 


GTDB  Event  ev120:  Steered  P 


GTDB  Event  ev124:  Steered  P 


GTDB  Event  ev125:  Steered  P  GTDB  Event  ev135:  Steered  P 


0  2  4  6  8  10  12  0  2  4  6  8  10  12 

Time  (sec)  Time  (sec) 


Figure  9.  Steered  Pn  waveforms  at  the  Sonseca  array  for  four  events  from  the  Puertollano  region. 
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Broadband  Pn/Lg 


o:  Earthquakes  x:  Explosions 


Figure  10.  Broadband  Pn/Lg  spectral  ratios  after  source  deconvolution  versus  derived  source  time 
function  for  events  at  NORESS,  Apatity,  Steigen,  Puertollano  Spain,  and  Hormuz. 
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Conclusions 


We  have  applied  a  source  deconvolution  process  to  the  Pn  and  Lg  waves  of  small  mining 
explosions  and  earthquakes  in  order  to  stabilize  the  spectral  ratio  estimates  and  improve 
discrimination  capability.  This  stabilized  spectral  ratio  is  then  used  in  conjunction  with  the  derived 
source  duration  for  event  identification.  This  procedure  has  been  tested  on  a  number  of  Ground 
Truth  Databases  which  include  events  from  Scandinavia,  central  Europe,  Spain,  and  the  Middle 
East.  We  have  found  that  the  scatter  in  the  spectral  ratios  and  their  parameterizations  is  reduced, 
leading  to  fewer  misidentified  events. 
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JAMES  E.  ZOLLWEG 
BOISE  STATE  UNIVERSITY 
GEOSCIENCES  DEPT. 

1910  UNIVERSITY  DRIVE 
BOISE,  ID  83725 


TERRY  WALLACE 
UNIVERSITY  OF  ARIZONA 
DEPARTMENT  OF  GEOSCIENCES 
BUILDING  #77 
TUCSON,  AZ  85721 

DANIEL  WEILL 
NSF 

EAR-785 

4201  WILSON  BLVD.,  ROOM  785 
ARLINGTON,  VA  22230 


JIAKANG  XIE 
COLUMBIA  UNIVERSITY 
LAMONT  DOHERTY  EARTH  OBSERV. 
ROUTE  9W 
PALISADES,  NY  10964 
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